#########################################################################
### Kline, Bankert, Levitan, Kraft. 2017.                             ###
### "Personality and Prosocial Behavior: A Multilevel Meta-Analysis"  ###
### - summary of results                                              ###
#########################################################################

## load packages etc
rm(list=ls())
library(rstan)
library(dplyr)
library(ggplot2)
library(stargazer)
library(lme4)

## set options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## load data and results from MLMA_bayes_estimation.R (wd is set in MASTER.R)
load("results.Rdata")

## labels for plots
lab_beta <- c("Extraversion","Openness","Agreeableness","Conscientousness","Neuroticism")
lab_gamma <- c("(Global) Intercept","Payment","Author: Hilbig","Author: Ben-Ner"
              ,"MBTI Personality Measure", "Cooperative")
lab_level <- c("Subject","Study")
lab_order <- rev(c(lab_beta, lab_gamma))
lab_study <- c("Brocklebank et al. 2011","Hilbig/Zettler 2009","Hilbig et al. 2012b"
               ,"Hilbig et al. 2012a","Schmitt et al. 2008","Swope et al. 2008"
               ,"Ben-Ner et al. 2004","Hirsh/Peterson 2009","Pothos et al. 2011"
               ,"Ben-Ner/Kramer 2010","Kurzban/Houser 2001","Koole et al. 2001")


### Summary statistics

stargazer(data.frame(dat1), title="Data for Model 1 (12 Studies)"
          , label = "tab:dat1", out="tab_dat1.tex", type="text")
stargazer(data.frame(dat2), title="Data for Model 2 (10 Studies)"
          , label = "tab:dat2", out="tab_dat2.tex", type="text")
stargazer(data.frame(dat3), title="Data for Model 3 (15 Studies)"
          , label = "tab:dat3", out="tab_dat3.tex", type="text")


### Original model results

df1 <- data.frame(summary(m1stan, par = c("beta", "gamma"))$summary) %>%
    mutate(Model = "Model 1", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta[1:4],lab_gamma)
         , Level = factor(var, labels = lab_level))
df2 <- data.frame(summary(m2stan, par = c("beta", "gamma"))$summary) %>%
    mutate(Model = "Model 2", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta,lab_gamma[-5])
         , Level = factor(var, labels = lab_level))
df3 <- data.frame(summary(m3stan, par = c("gamma"))$summary) %>%
    mutate(Model = "Model 3", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = lab_gamma[c(1,2,6)]
         , Level = factor(lab_level[2], levels = lab_level))
df_plot <- bind_rows(df1, df2, df3) %>% mutate(Variable = factor(Variable, levels = lab_order))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5., shape = Level)) +
    coord_flip() + geom_pointrange() + theme_bw() + facet_wrap(~Model) + ylab("Estimate") +
    geom_hline(yintercept = 0, col = "grey") + theme(legend.position = "bottom")
ggsave("fig1_original.pdf", width = 7, height = 4)


### Robustness check: Rescaling pro-social behavior (model 1-3)

df1 <- data.frame(summary(m1svo2stan, par = c("beta", "gamma"))$summary) %>%
  mutate(Model = "Model 1", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta[1:4],lab_gamma)
         , Level = factor(var, labels = lab_level))
df2 <- data.frame(summary(m2svo2stan, par = c("beta", "gamma"))$summary) %>%
  mutate(Model = "Model 2", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta,lab_gamma[-5])
         , Level = factor(var, labels = lab_level))
df3 <- data.frame(summary(m3svo2stan, par = c("gamma"))$summary) %>%
  mutate(Model = "Model 3", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = lab_gamma[c(1,2,6)]
         , Level = factor(lab_level[2], levels = lab_level))
df_plot <- bind_rows(df1, df2, df3) %>% mutate(Variable = factor(Variable, levels = lab_order))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5., shape = Level)) +
  coord_flip() + geom_pointrange() + theme_bw() + facet_wrap(~Model) + ylab("Estimate") +
  geom_hline(yintercept = 0, col = "grey") + theme(legend.position = "bottom")
ggsave("fig2_svo.pdf", width = 7, height = 4)


### Robustness check: Rescaling pro-social behavior / SVO variable to dummy (model 1-3)

df1 <- data.frame(summary(m1svoD1stan, par = c("beta", "gamma"))$summary) %>%
  rbind(data.frame(summary(m1svoD2stan, par = c("beta", "gamma"))$summary)) %>%
  mutate(Model = "Model 1", var = gsub("\\[\\d*\\]\\d*","",rownames(.))
         , Variable = rep(c(lab_beta[1:4],lab_gamma), 2)
         , Level = factor(var, labels = lab_level)
         , svo = rep(c("Dichotomized DV (svoD1)", "Dichotomized DV (svoD2)"), each=nrow(.)/2))
df2 <- data.frame(summary(m2svoD1stan, par = c("beta", "gamma"))$summary) %>%
  rbind(data.frame(summary(m2svoD2stan, par = c("beta", "gamma"))$summary)) %>%
  mutate(Model = "Model 2", var = gsub("\\[\\d*\\]\\d*","",rownames(.))
         , Variable = rep(c(lab_beta,lab_gamma[-5]), 2)
         , Level = factor(var, labels = lab_level)
         , svo = rep(c("Dichotomized DV (svoD1)", "Dichotomized DV (svoD2)"), each=nrow(.)/2))
df3 <- data.frame(summary(m3svoD1stan, par = "gamma")$summary) %>%
  rbind(data.frame(summary(m3svoD2stan, par = "gamma")$summary)) %>%
  mutate(Model = "Model 3", var = gsub("\\[\\d*\\]\\d*","",rownames(.))
         , Variable = rep(lab_gamma[c(1,2,6)], 2)
         , Level = factor(lab_level[2], levels = lab_level)
         , svo = rep(c("Dichotomized DV (svoD1)", "Dichotomized DV (svoD2)"), each=nrow(.)/2))
df_plot <- bind_rows(df1, df2, df3) %>% mutate(Variable = factor(Variable, levels = lab_order))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5., shape = Level)) +
  coord_flip() + geom_pointrange() + theme_bw() + facet_grid(svo~Model) + ylab("Estimate") +
  geom_hline(yintercept = 0, col = "grey") + theme(legend.position = "bottom")
ggsave("fig3_logit.pdf", width = 7, height = 5)


### Additional robustness checks (focusing on model 1)

## 1) varying study inclusion
df_prep <- function(x, i = NULL){
    out <- data.frame(summary(x, par = c("beta", "gamma"))$summary) %>%
        mutate(var = gsub("\\[\\d*\\]","",rownames(.))
             , Variable = c(lab_beta[1:4],lab_gamma)
             , Level = factor(var, labels = lab_level))
    return(out)
}
df_plot <- m4stan %>% lapply(df_prep) %>% bind_rows() %>%
    mutate(Study = rep(lab_study, each = 10)
         , Variable = factor(Variable, levels = lab_order))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5., shape = Level)) +
    coord_flip() + geom_pointrange() + theme_bw() + facet_wrap(~Study, ncol = 3) + ylab("Estimate") +
    geom_hline(yintercept = 0, col = "grey") + theme(legend.position = "bottom")
ggsave("fig4_jackknife.pdf", width = 7, height = 8)

## 2) random slope models
df_plot <- data.frame(summary(m5stan, par = c("beta_extro","beta_open","beta_agree"
                                             ,"beta_conscien", "mu_beta","gamma"))$summary) %>%
    mutate(var = rownames(.)
         , Variable = factor(c(rep(lab_beta[1:4], each = 12),lab_beta[1:4],lab_gamma)
                           , levels = rev(lab_order))
         , Study = factor(c(rep(lab_study, 4), rep("Overall", 4), rep("Study-level", 6))
                        , levels = c(sort(lab_study), "Overall", "Study-level")))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5.)) +
    geom_pointrange() + theme_bw() + facet_wrap(~Study, scales="free_x") +
    ylab("Estimate") + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    geom_hline(yintercept = 0, col = "grey")
ggsave("fig5_slopes.pdf", width = 9, height = 11)


### Robustness check: drop cooperation dummy (model 1-3)

df1 <- data.frame(summary(m1coop0stan, par = c("beta", "gamma"))$summary) %>%
  mutate(Model = "Model 1", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta[1:4],lab_gamma[-6])
         , Level = factor(var, labels = lab_level))
df2 <- data.frame(summary(m2coop0stan, par = c("beta", "gamma"))$summary) %>%
  mutate(Model = "Model 2", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = c(lab_beta,lab_gamma[-5:-6])
         , Level = factor(var, labels = lab_level))
df3 <- data.frame(summary(m3coop0stan, par = c("gamma"))$summary) %>%
  mutate(Model = "Model 3", var = gsub("\\[\\d*\\]","",rownames(.))
         , Variable = lab_gamma[c(1,2)]
         , Level = factor(lab_level[2], levels = lab_level))
df_plot <- bind_rows(df1, df2, df3) %>% mutate(Variable = factor(Variable, levels = lab_order))

ggplot(df_plot, aes(x = Variable, y = mean, ymin = X2.5., ymax = X97.5., shape = Level)) +
  coord_flip() + geom_pointrange() + theme_bw() + facet_wrap(~Model) + ylab("Estimate") +
  geom_hline(yintercept = 0, col = "grey") + theme(legend.position = "bottom")
ggsave("fig6_cooperative.pdf", width = 7, height = 4)


### Frequentist model estimation

## Model 1: Big 5 + MBTI (12 studies)
m1 <- lmer(svo1 ~ extro + open + agree + conscien +
             payment + hilbig + benner + mbti + coop +
             (1 | study), data = dat1)

## Model 2: Big 5 (10 studies)
m2 <- lmer(svo1 ~ extro + open + agree + conscien + neuro +
             payment + hilbig + benner + coop +
             (1 | study), data = dat2)

## Model 3: Incentivization (15 studies)
m3 <- lmer(svo1 ~ payment + coop +
             (1 | study), data = dat3)

## model output
stargazer(m1, m2, m3, ci = T
          , covariate.labels = c(lab_beta, lab_gamma[-1], lab_gamma), label = "tab:mfreq"
          , title = "Frequentist replication of main models (multilevel regression)"
          , dep.var.labels = "Prosociality", keep.stat = "n"
          , star.cutoffs = NA,omit.table.layout = "n", out="tab_mfreq.tex", type="text")


### Traceplots for original models (in appendix)
# NOTE: traceplot via rstan produce ggplot error on some platforms
# (check package updates if that is the case!)

# produce traceplots
rstan::traceplot(m1stan, ncol = 3) + ggtitle("Model 1")
try(ggsave("fig_trace1.pdf",height = 6, width=9), silent = T)
try(ggsave("fig_trace1.png",height = 6, width=9), silent = T)
rstan::traceplot(m2stan, ncol = 3) + ggtitle("Model 2")
try(ggsave("fig_trace2.pdf",height = 6, width=9), silent = T)
try(ggsave("fig_trace2.png",height = 6, width=9), silent = T)
rstan::traceplot(m3stan, pars = "gamma", ncol = 3) + ggtitle("Model 3")
try(ggsave("fig_trace3.pdf",height = 2, width=9), silent = T)
try(ggsave("fig_trace3.png",height = 2, width=9), silent = T)


### Traceplots via coda (in case of ggplot error)

# install coda
if(!"coda" %in% installed.packages()) {
  install.packages("coda")
}

# load coda
library(coda)

# produce traceplots using coda
pdf("fig_trace1_coda.pdf",height = 9, width=9)
par(mfrow=c(4,3))
coda::traceplot(As.mcmc.list(m1stan, pars=c("beta","gamma")))
dev.off()

pdf("fig_trace2_coda.pdf",height = 9, width=9)
par(mfrow=c(4,3))
coda::traceplot(As.mcmc.list(m2stan, pars=c("beta","gamma")))
dev.off()

pdf("fig_trace3_coda.pdf",height = 3, width=9)
par(mfrow=c(1,3))
coda::traceplot(As.mcmc.list(m3stan, pars="gamma"))
dev.off()
